
/*
 *
 * HexMetric.cpp contains quality calculations for hexes
 *
 * Copyright (C) 2003 Sandia National Laboratories <cubit@sandia.gov>
 *
 * This file is part of VERDICT
 *
 * This copy of VERDICT is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * VERDICT is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */



#include "VerdictVector.h"
#include "V_GaussIntegration.h"
#include "verdict_defines.h"
#include <memory.h>
#include <algorithm>

//! the average volume of a hex
VERDICT_REAL verdict_hex_size = 0;

//! weights based on the average size of a hex
int v_hex_get_weight( VerdictVector &v1, 
    VerdictVector &v2,
    VerdictVector &v3 )
{

  if( verdict_hex_size == 0)
    return 0;

  v1.set(1,0,0);
  v2.set(0,1,0);
  v3.set(0,0,1);

  double scale = pow(verdict_hex_size/ (v1 % (v2 * v3 )), 0.33333333333); 
  v1 *= scale;
  v2 *= scale;
  v3 *= scale;

  return 1;
}

//! returns the average volume of a hex
C_FUNC_DEF void v_set_hex_size( VERDICT_REAL size )
{
  verdict_hex_size = size;
}

#define make_hex_nodes(coord, pos)                      \
  for (int mhcii = 0; mhcii < 8; mhcii++ )              \
  {                                                     \
    pos[mhcii].set( coord[mhcii][0],                    \
        coord[mhcii][1],                                \
        coord[mhcii][2] );                              \
  }


#define make_edge_length_squares(edges, lengths)        \
{                                                       \
  for(int melii=0; melii<12; melii++)                   \
    lengths[melii] = edges[melii].length_squared();     \
}
  

//! make VerdictVectors from coordinates
void make_hex_edges( VERDICT_REAL coordinates[][3], VerdictVector edges[12] )
{
  edges[0].set(
      coordinates[1][0] - coordinates[0][0],
      coordinates[1][1] - coordinates[0][1],
      coordinates[1][2] - coordinates[0][2]
      );
  edges[1].set(
      coordinates[2][0] - coordinates[1][0],
      coordinates[2][1] - coordinates[1][1],
      coordinates[2][2] - coordinates[1][2]
      );
  edges[2].set(
      coordinates[3][0] - coordinates[2][0],
      coordinates[3][1] - coordinates[2][1],
      coordinates[3][2] - coordinates[2][2]
      );
  edges[3].set(
      coordinates[0][0] - coordinates[3][0],
      coordinates[0][1] - coordinates[3][1],
      coordinates[0][2] - coordinates[3][2]
      );
  edges[4].set(
      coordinates[5][0] - coordinates[4][0],
      coordinates[5][1] - coordinates[4][1],
      coordinates[5][2] - coordinates[4][2]
      );
  edges[5].set(
      coordinates[6][0] - coordinates[5][0],
      coordinates[6][1] - coordinates[5][1],
      coordinates[6][2] - coordinates[5][2]
      );
  edges[6].set(
      coordinates[7][0] - coordinates[6][0],
      coordinates[7][1] - coordinates[6][1],
      coordinates[7][2] - coordinates[6][2]
      );
  edges[7].set(
      coordinates[4][0] - coordinates[7][0],
      coordinates[4][1] - coordinates[7][1],
      coordinates[4][2] - coordinates[7][2]
      );
  edges[8].set(
      coordinates[4][0] - coordinates[0][0],
      coordinates[4][1] - coordinates[0][1],
      coordinates[4][2] - coordinates[0][2]
      );
  edges[9].set(
      coordinates[5][0] - coordinates[1][0],
      coordinates[5][1] - coordinates[1][1],
      coordinates[5][2] - coordinates[1][2]
      );
  edges[10].set(
      coordinates[6][0] - coordinates[2][0],
      coordinates[6][1] - coordinates[2][1],
      coordinates[6][2] - coordinates[2][2]
      );
  edges[11].set(
      coordinates[7][0] - coordinates[3][0],
      coordinates[7][1] - coordinates[3][1],
      coordinates[7][2] - coordinates[3][2]
      );
}

/*!
  localizes hex coordinates
*/
void localize_hex_coordinates(VERDICT_REAL coordinates[][3], VerdictVector position[8] )
{

  int ii;
  for ( ii = 0; ii < 8; ii++ )
  {
    position[ii].set( coordinates[ii][0],
        coordinates[ii][1],
        coordinates[ii][2] );
  }
  
  // ... Make centroid of element the center of coordinate system
  VerdictVector point_1256 = position[1];
  point_1256 += position[2];
  point_1256 += position[5];
  point_1256 += position[6];

  VerdictVector point_0374 = position[0];
  point_0374 += position[3];
  point_0374 += position[7];
  point_0374 += position[4];

  VerdictVector centroid = point_1256;
  centroid += point_0374;
  centroid /= 8.0;

  int i;
  for ( i = 0; i < 8; i++)
    position[i] -= centroid;

  // ... Rotate element such that center of side 1-2 and 0-3 define X axis
  double DX = point_1256.x() - point_0374.x();
  double DY = point_1256.y() - point_0374.y();
  double DZ = point_1256.z() - point_0374.z();

  double AMAGX = sqrt(DX*DX + DZ*DZ);
  double AMAGY = sqrt(DX*DX + DY*DY + DZ*DZ);
  double FMAGX = AMAGX == 0.0 ? 1.0 : 0.0;
  double FMAGY = AMAGY == 0.0 ? 1.0 : 0.0;

  double CZ = DX / (AMAGX + FMAGX) + FMAGX;
  double SZ = DZ / (AMAGX + FMAGX);
  double CY = AMAGX / (AMAGY + FMAGY) + FMAGY;
  double SY = DY / (AMAGY + FMAGY);

  double temp;
 
  for (i = 0; i < 8; i++) 
  {
    temp =  CY * CZ * position[i].x() + CY * SZ * position[i].z() +
      SY * position[i].y();
    position[i].y( -SY * CZ * position[i].x() - SY * SZ * position[i].z() +
        CY * position[i].y());
    position[i].z( -SZ * position[i].x() + CZ * position[i].z());
    position[i].x(temp);
  }


  // ... Now, rotate about Y
  VerdictVector delta = -position[0];
  delta -= position[1];
  delta += position[2];
  delta += position[3];
  delta -= position[4];
  delta -= position[5];
  delta += position[6];
  delta += position[7];

  DX = delta.x();
  DY = delta.y();
  DZ = delta.z();

  AMAGY = sqrt(DY*DY + DZ*DZ);
  FMAGY = AMAGY == 0.0 ? 1.0 : 0.0;

  double CX = DY / (AMAGY + FMAGY) + FMAGY;
  double SX = DZ / (AMAGY + FMAGY);
  
  for (i = 0; i < 8; i++) 
  {
    temp =  CX * position[i].y() + SX * position[i].z();
    position[i].z(-SX * position[i].y() + CX * position[i].z());
    position[i].y(temp);
  }
}


double safe_ratio3( const double numerator, 
    const double denominator,
    const double max_ratio )
{
    // this filter is essential for good running time in practice
  double return_value = 0.0;

  const double filter_n = max_ratio * 1.0e-16;
  const double filter_d = 1.0e-16;
  if ( fabs( numerator ) <= filter_n &&
      fabs( denominator ) >= filter_d )
    return_value = numerator / denominator;
  
  else
    return_value = fabs(numerator) / max_ratio >= fabs(denominator) ?
      ( (numerator >= 0.0 && denominator >= 0.0) ||
        (numerator < 0.0 && denominator < 0.0) ?
        max_ratio : -max_ratio )
      : numerator / denominator;
  
  return return_value;
}


double safe_ratio( const double numerator, 
    const double denominator )
{

  double return_value = 0.0;
  const double filter_n = VERDICT_DBL_MAX; 
  const double filter_d = VERDICT_DBL_MIN;
  if ( fabs( numerator ) <= filter_n &&
      fabs( denominator ) >= filter_d )
    return_value = numerator / denominator;
  
  else
    return_value = VERDICT_DBL_MAX;

  return return_value;
  
}



double condition_comp( const VerdictVector &xxi, 
    const VerdictVector &xet, 
    const VerdictVector &xze )
{
  double det =  xxi % (xet * xze);
  
  if( det <= VERDICT_DBL_MIN) { return VERDICT_DBL_MAX; }
  
  
  double term1 = xxi % xxi + xet % xet + xze % xze;
  double term2 = (( xxi*xet ) % ( xxi*xet )) + (( xet*xze ) % ( xet*xze )) + (( xze*xxi ) % ( xze*xxi ));
  
  return sqrt( term1 * term2 ) / det;
}



double oddy_comp( const VerdictVector &xxi, 
    const VerdictVector &xet, 
    const VerdictVector &xze )
{
  static const double third=1.0/3.0;
  
  double g11, g12, g13, g22, g23, g33, rt_g;
  
  g11 = xxi % xxi;
  g12 = xxi % xet;
  g13 = xxi % xze;
  g22 = xet % xet;
  g23 = xet % xze;
  g33 = xze % xze;
  rt_g = xxi % (xet * xze);
  
  double oddy_metric;
  if( rt_g > VERDICT_DBL_MIN ) 
  {
    double norm_G_squared = g11*g11 + 2.0*g12*g12 + 2.0*g13*g13 + g22*g22 + 2.0*g23*g23 +g33*g33;
    
    double norm_J_squared = g11 + g22 + g33;
    
    oddy_metric = ( norm_G_squared - third*norm_J_squared*norm_J_squared ) / pow( rt_g, 4.*third );
    
  }
  
  else
    oddy_metric = VERDICT_DBL_MAX;
  
  return oddy_metric;
  
}


//! calcualates edge lengths of a hex
VERDICT_REAL hex_edge_length(int max_min, VERDICT_REAL coordinates[][3])
{
  double temp[3], edge[12];
  int i;
  
  //lengths^2 of edges
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[1][i] - coordinates[0][i];
    temp[i] = temp[i] * temp[i];
  }
  edge[0] = sqrt( temp[0] + temp[1] + temp[2] );
  
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[2][i] - coordinates[1][i];
    temp[i] = temp[i] * temp[i];
  }
  edge[1] = sqrt( temp[0] + temp[1] + temp[2] );
  
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[3][i] - coordinates[2][i];
    temp[i] = temp[i] * temp[i];
  }
  edge[2] = sqrt( temp[0] + temp[1] + temp[2] );
  
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[0][i] - coordinates[3][i];
    temp[i] = temp[i] * temp[i];
  }
  edge[3] = sqrt( temp[0] + temp[1] + temp[2] );
  
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[5][i] - coordinates[4][i];
    temp[i] = temp[i] * temp[i];
  }
  edge[4] = sqrt( temp[0] + temp[1] + temp[2] );
  
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[6][i] - coordinates[5][i];
    temp[i] = temp[i] * temp[i];
  }
  edge[5] = sqrt( temp[0] + temp[1] + temp[2] );
  
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[7][i] - coordinates[6][i];
    temp[i] = temp[i] * temp[i];
  }
  edge[6] = sqrt( temp[0] + temp[1] + temp[2] );
  
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[4][i] - coordinates[7][i];
    temp[i] = temp[i] * temp[i];
  }
  edge[7] = sqrt( temp[0] + temp[1] + temp[2] );
  
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[4][i] - coordinates[0][i];
    temp[i] = temp[i] * temp[i];
  }
  edge[8] = sqrt( temp[0] + temp[1] + temp[2] );
  
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[5][i] - coordinates[1][i];
    temp[i] = temp[i] * temp[i];
  }
  edge[9] = sqrt( temp[0] + temp[1] + temp[2] );
  
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[6][i] - coordinates[2][i];
    temp[i] = temp[i] * temp[i];
  }
  edge[10] = sqrt( temp[0] + temp[1] + temp[2] );
  
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[7][i] - coordinates[3][i];
    temp[i] = temp[i] * temp[i];
  }
  edge[11] = sqrt( temp[0] + temp[1] + temp[2] );
  
  double _edge = edge[0];
  
  if(max_min == 0)
  {
    for( i = 1; i<12; i++) 
      _edge = std::min( _edge, edge[i] ); 
    return (VERDICT_REAL)_edge;
  }  
  else
  {
    for( i = 1; i<12; i++) 
      _edge = std::max( _edge, edge[i] );
    return (VERDICT_REAL)_edge;
  }
  
}


VERDICT_REAL diag_length(int max_min, VERDICT_REAL coordinates[][3]) 
{
  double temp[3], diag[4];
  int i = 0;
  
  //lengths^2  f diag nals
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[6][i] - coordinates[0][i];
    temp[i] = temp[i] * temp[i];
  }
  diag[0] = sqrt( temp[0] + temp[1] + temp[2] );

  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[4][i] - coordinates[2][i];
    temp[i] = temp[i] * temp[i];
  }
  diag[1] = sqrt( temp[0] + temp[1] + temp[2] );

  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[7][i] - coordinates[1][i];
    temp[i] = temp[i] * temp[i];
  }
  diag[2] = sqrt( temp[0] + temp[1] + temp[2] );

  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[5][i] - coordinates[3][i];
    temp[i] = temp[i] * temp[i];
  }
  diag[3] = sqrt( temp[0] + temp[1] + temp[2] );

  double diagonal = diag[0];
  if(max_min == 0 )  //Return min diagonal
  { 
    for( i = 1; i<4; i++)
      diagonal = std::min( diagonal, diag[i] );
    return (VERDICT_REAL)diagonal;
  }
  else          //Return max diagonal
  {
    for( i = 1; i<4; i++)
      diagonal = std::max( diagonal, diag[i] );
    return (VERDICT_REAL)diagonal;  
  }

}

//! calculates efg values
VerdictVector calc_hex_efg( int efg_index, VerdictVector coordinates[8])
{
  
  VerdictVector efg;
  
  switch(efg_index) {
    
    case 1:
      efg =  coordinates[1];
      efg += coordinates[2];
      efg += coordinates[5];
      efg += coordinates[6];
      efg -= coordinates[0];
      efg -= coordinates[3];
      efg -= coordinates[4];
      efg -= coordinates[7];
      break;
      
    case 2:
      efg =  coordinates[2];
      efg += coordinates[3];
      efg += coordinates[6];
      efg += coordinates[7];
      efg -= coordinates[0];
      efg -= coordinates[1];
      efg -= coordinates[4];
      efg -= coordinates[5];
      break;
      
    case 3:
      efg =  coordinates[4];
      efg += coordinates[5];
      efg += coordinates[6];
      efg += coordinates[7];
      efg -= coordinates[0];
      efg -= coordinates[1];
      efg -= coordinates[2];
      efg -= coordinates[3];
      break;

    case 12:
      efg =  coordinates[0];
      efg += coordinates[2];
      efg += coordinates[4];
      efg += coordinates[6];
      efg -= coordinates[1];
      efg -= coordinates[3];
      efg -= coordinates[5];
      efg -= coordinates[7];
      break;
      
    case 13:
      efg =  coordinates[0];
      efg += coordinates[3];
      efg += coordinates[5];
      efg += coordinates[6];
      efg -= coordinates[1];
      efg -= coordinates[2];
      efg -= coordinates[4];
      efg -= coordinates[7];
      break;

   case 23:
      efg =  coordinates[0];
      efg += coordinates[1];
      efg += coordinates[6];
      efg += coordinates[7];
      efg -= coordinates[2];
      efg -= coordinates[3];
      efg -= coordinates[4];
      efg -= coordinates[5];
      break;

   case 123:
      efg =  coordinates[0];
      efg += coordinates[2];
      efg += coordinates[5];
      efg += coordinates[7];
      efg -= coordinates[1];
      efg -= coordinates[5];
      efg -= coordinates[6];
      efg -= coordinates[2];
      break;
      
   default:
      efg.set(0,0,0);
      
  }
  
  return efg;
}

/*!
  the aspect ratio of a hexahedral

  Maximum edge length ratios at hex center
*/
C_FUNC_DEF VERDICT_REAL v_hex_aspect (int /*num_nodes*/, VERDICT_REAL coordinates[][3])
{
  double aspect;
  VerdictVector node_pos[8];
  make_hex_nodes ( coordinates, node_pos );
  
  double aspect_12, aspect_13, aspect_23;
  
  VerdictVector efg1 = calc_hex_efg( 1, node_pos);
  VerdictVector efg2 = calc_hex_efg( 2, node_pos);
  VerdictVector efg3 = calc_hex_efg( 3, node_pos);

  double mag_efg1 = efg1.length();
  double mag_efg2 = efg2.length();
  double mag_efg3 = efg3.length();
  
  aspect_12 = safe_ratio( std::max( mag_efg1, mag_efg2 ) , std::min( mag_efg1, mag_efg2 ) );
  aspect_13 = safe_ratio( std::max( mag_efg1, mag_efg3 ) , std::min( mag_efg1, mag_efg3 ) );
  aspect_23 = safe_ratio( std::max( mag_efg2, mag_efg3 ) , std::min( mag_efg2, mag_efg3 ) );

  aspect = std::max( aspect_12, std::max( aspect_13, aspect_23 ) );

  if( aspect > 0 )
    return (VERDICT_REAL) std::min( aspect, VERDICT_DBL_MAX );
  return (VERDICT_REAL) std::max( aspect, -VERDICT_DBL_MAX );
 
}

/*!
  skew of a hex

  Maximum ||cosA|| where A is the angle between edges at hex center.
*/
C_FUNC_DEF VERDICT_REAL v_hex_skew( int /*num_nodes*/, VERDICT_REAL coordinates[][3] )
{
  VerdictVector node_pos[8];
  make_hex_nodes ( coordinates, node_pos );
  
  double skew_1, skew_2, skew_3;
  
  VerdictVector efg1 = calc_hex_efg( 1, node_pos);
  VerdictVector efg2 = calc_hex_efg( 2, node_pos);
  VerdictVector efg3 = calc_hex_efg( 3, node_pos);
 
  if( efg1.normalize() <= VERDICT_DBL_MIN )
    return VERDICT_DBL_MAX;
  if( efg2.normalize() <= VERDICT_DBL_MIN )
    return VERDICT_DBL_MAX;
  if( efg3.normalize() <= VERDICT_DBL_MIN )
    return VERDICT_DBL_MAX;

  skew_1 = fabs(efg1 % efg2);
  skew_2 = fabs(efg1 % efg3);
  skew_3 = fabs(efg2 % efg3);

  double skew = (std::max( skew_1, std::max( skew_2, skew_3 ) ));

  if( skew > 0 )
    return (VERDICT_REAL) std::min( skew, VERDICT_DBL_MAX );
  return (VERDICT_REAL) std::max( skew, -VERDICT_DBL_MAX );
}

/*!
  taper of a hex

  Maximum ratio of lengths derived from opposite edges.
*/
C_FUNC_DEF VERDICT_REAL v_hex_taper( int /*num_nodes*/, VERDICT_REAL coordinates[][3] )
{
  VerdictVector node_pos[8];
  make_hex_nodes ( coordinates, node_pos );
  
  VerdictVector efg1 = calc_hex_efg( 1, node_pos);
  VerdictVector efg2 = calc_hex_efg( 2, node_pos);
  VerdictVector efg3 = calc_hex_efg( 3, node_pos);

  VerdictVector efg12 = calc_hex_efg( 12, node_pos);
  VerdictVector efg13 = calc_hex_efg( 13, node_pos);
  VerdictVector efg23 = calc_hex_efg( 23, node_pos);

  double taper_1 = fabs( safe_ratio( efg12.length(), std::min( efg1.length(), efg2.length())));
  double taper_2 = fabs( safe_ratio( efg13.length(), std::min( efg1.length(), efg3.length())));
  double taper_3 = fabs( safe_ratio( efg23.length(), std::min( efg2.length(), efg3.length())));

  double taper = (VERDICT_REAL)std::max(taper_1, std::max(taper_2, taper_3));  

  if( taper > 0 )
    return (VERDICT_REAL) std::min( taper, VERDICT_DBL_MAX );
  return (VERDICT_REAL) std::max( taper, -VERDICT_DBL_MAX );
}

/*!
  volume of a hex

  Jacobian at hex center
*/
C_FUNC_DEF VERDICT_REAL v_hex_volume( int /*num_nodes*/, VERDICT_REAL coordinates[][3] )
{
  VerdictVector node_pos[8];
  make_hex_nodes ( coordinates, node_pos );
  
  VerdictVector efg1 = calc_hex_efg( 1, node_pos);
  VerdictVector efg2 = calc_hex_efg( 2, node_pos);
  VerdictVector efg3 = calc_hex_efg( 3, node_pos);

  double volume;
  volume = (VERDICT_REAL) (efg1 % (efg2 * efg3))/64.0;

  if( volume > 0 )
    return (VERDICT_REAL) std::min( volume, VERDICT_DBL_MAX );
  return (VERDICT_REAL) std::max( volume, -VERDICT_DBL_MAX );
}

/*!
  stretch of a hex

  sqrt(3) * minimum edge length / maximum diagonal length
*/
C_FUNC_DEF VERDICT_REAL v_hex_stretch( int /*num_nodes*/, VERDICT_REAL coordinates[][3] )
{
  static const double HEX_STRETCH_SCALE_FACTOR = sqrt(3.0);
  
  double min_edge = hex_edge_length( 0, coordinates );
  double max_diag = diag_length( 1, coordinates );  
  
  double stretch = HEX_STRETCH_SCALE_FACTOR * safe_ratio(min_edge, max_diag);

  if( stretch > 0 )
    return (VERDICT_REAL) std::min( stretch, VERDICT_DBL_MAX );
  return (VERDICT_REAL) std::max( stretch, -VERDICT_DBL_MAX );
}

/*!
  Min diagonal length of a hex
*/
C_FUNC_DEF VERDICT_REAL v_hex_min_diagonal( int /*num_nodes*/, VERDICT_REAL coordinates[][3] )
{
  double min_diag = diag_length( 0, coordinates );
  
  return (VERDICT_REAL)min_diag;
}

/*!
  Max diagonal length of a hex
*/
C_FUNC_DEF VERDICT_REAL v_hex_max_diagonal( int /*num_nodes*/, VERDICT_REAL coordinates[][3] )
{
  double max_diag = diag_length( 1, coordinates );
  
  return (VERDICT_REAL)max_diag;
}


/*!
  diagonal ratio of a hex
  
  Minimum diagonal length / maximum diagonal length
*/
C_FUNC_DEF VERDICT_REAL v_hex_diagonal_ratio( int /*num_nodes*/, VERDICT_REAL coordinates[][3] )
{
  
  double min_diag = diag_length( 0, coordinates ); 
  double max_diag = diag_length( 1, coordinates );
  
  double diagonal = safe_ratio( min_diag, max_diag);

  if( diagonal > 0 )
    return (VERDICT_REAL) std::min( diagonal, VERDICT_DBL_MAX );
  return (VERDICT_REAL) std::max( diagonal, -VERDICT_DBL_MAX );
}

#define SQR(x) ((x) * (x))

/*!
  dimension of a hex

  Pronto-specific characteristic length for stable time step calculation. 
  Char_length = Volume / 2 grad Volume
*/
C_FUNC_DEF VERDICT_REAL v_hex_dimension( int /*num_nodes*/, VERDICT_REAL coordinates[][3] )
{
  
  double gradop[9][4];

  double x1 = coordinates[0][0];
  double x2 = coordinates[1][0];
  double x3 = coordinates[2][0];
  double x4 = coordinates[3][0];
  double x5 = coordinates[4][0];
  double x6 = coordinates[5][0];
  double x7 = coordinates[6][0];
  double x8 = coordinates[7][0];

  double y1 = coordinates[0][1];
  double y2 = coordinates[1][1];
  double y3 = coordinates[2][1];
  double y4 = coordinates[3][1];
  double y5 = coordinates[4][1];
  double y6 = coordinates[5][1];
  double y7 = coordinates[6][1];
  double y8 = coordinates[7][1];

  double z1 = coordinates[0][2];
  double z2 = coordinates[1][2];
  double z3 = coordinates[2][2];
  double z4 = coordinates[3][2];
  double z5 = coordinates[4][2];
  double z6 = coordinates[5][2];
  double z7 = coordinates[6][2];
  double z8 = coordinates[7][2];

  double z24 = z2 - z4;
  double z52 = z5 - z2;
  double z45 = z4 - z5;
  gradop[1][1] = ( y2*(z6-z3-z45) + y3*z24 + y4*(z3-z8-z52)
      + y5*(z8-z6-z24) + y6*z52 + y8*z45 ) / 12.0;
         
  double z31 = z3 - z1;
  double z63 = z6 - z3;
  double z16 = z1 - z6;
  gradop[2][1] = ( y3*(z7-z4-z16) + y4*z31 + y1*(z4-z5-z63)
      + y6*(z5-z7-z31) + y7*z63 + y5*z16 ) / 12.0;

  double z42 = z4 - z2;
  double z74 = z7 - z4;
  double z27 = z2 - z7;
  gradop[3][1] = ( y4*(z8-z1-z27) + y1*z42 + y2*(z1-z6-z74)
      + y7*(z6-z8-z42) + y8*z74 + y6*z27 ) / 12.0;
             
  double z13 = z1 - z3;
  double z81 = z8 - z1;
  double z38 = z3 - z8;
  gradop[4][1] = ( y1*(z5-z2-z38) + y2*z13 + y3*(z2-z7-z81)
      + y8*(z7-z5-z13) + y5*z81 + y7*z38 ) / 12.0;
   
  double z86 = z8 - z6;
  double z18 = z1 - z8;
  double z61 = z6 - z1;
  gradop[5][1] = ( y8*(z4-z7-z61) + y7*z86 + y6*(z7-z2-z18)
      + y1*(z2-z4-z86) + y4*z18 + y2*z61 ) / 12.0;
     
  double z57 = z5 - z7;
  double z25 = z2 - z5;
  double z72 = z7 - z2;
  gradop[6][1] = ( y5*(z1-z8-z72) + y8*z57 + y7*(z8-z3-z25)
      + y2*(z3-z1-z57) + y1*z25 + y3*z72 ) / 12.0;
       
  double z68 = z6 - z8;
  double z36 = z3 - z6;
  double z83 = z8 - z3;
  gradop[7][1] = ( y6*(z2-z5-z83) + y5*z68 + y8*(z5-z4-z36)
      + y3*(z4-z2-z68) + y2*z36 + y4*z83 ) / 12.0;
         
  double z75 = z7 - z5;
  double z47 = z4 - z7;
  double z54 = z5 - z4;
  gradop[8][1] = ( y7*(z3-z6-z54) + y6*z75 + y5*(z6-z1-z47)
      + y4*(z1-z3-z75) + y3*z47 + y1*z54 ) / 12.0;
           
  double x24 = x2 - x4;
  double x52 = x5 - x2;
  double x45 = x4 - x5;
  gradop[1][2] = ( z2*(x6-x3-x45) + z3*x24 + z4*(x3-x8-x52)
      + z5*(x8-x6-x24) + z6*x52 + z8*x45 ) / 12.0;

  double x31 = x3 - x1;
  double x63 = x6 - x3;
  double x16 = x1 - x6;
  gradop[2][2] = ( z3*(x7-x4-x16) + z4*x31 + z1*(x4-x5-x63)
      + z6*(x5-x7-x31) + z7*x63 + z5*x16 ) / 12.0;

  double x42 = x4 - x2;
  double x74 = x7 - x4;
  double x27 = x2 - x7;
  gradop[3][2] = ( z4*(x8-x1-x27) + z1*x42 + z2*(x1-x6-x74)
      + z7*(x6-x8-x42) + z8*x74 + z6*x27 ) / 12.0;

  double x13 = x1 - x3;
  double x81 = x8 - x1;
  double x38 = x3 - x8;
  gradop[4][2] = ( z1*(x5-x2-x38) + z2*x13 + z3*(x2-x7-x81)
      + z8*(x7-x5-x13) + z5*x81 + z7*x38 ) / 12.0;

  double x86 = x8 - x6;
  double x18 = x1 - x8;
  double x61 = x6 - x1;
  gradop[5][2] = ( z8*(x4-x7-x61) + z7*x86 + z6*(x7-x2-x18)
      + z1*(x2-x4-x86) + z4*x18 + z2*x61 ) / 12.0;
                     
  double x57 = x5 - x7;
  double x25 = x2 - x5;
  double x72 = x7 - x2;
  gradop[6][2] = ( z5*(x1-x8-x72) + z8*x57 + z7*(x8-x3-x25)
      + z2*(x3-x1-x57) + z1*x25 + z3*x72 ) / 12.0;

  double x68 = x6 - x8;
  double x36 = x3 - x6;
  double x83 = x8 - x3;
  gradop[7][2] = ( z6*(x2-x5-x83) + z5*x68 + z8*(x5-x4-x36)
      + z3*(x4-x2-x68) + z2*x36 + z4*x83 ) / 12.0;

  double x75 = x7 - x5;
  double x47 = x4 - x7;
  double x54 = x5 - x4;
  gradop[8][2] = ( z7*(x3-x6-x54) + z6*x75 + z5*(x6-x1-x47)
      + z4*(x1-x3-x75) + z3*x47 + z1*x54 ) / 12.0;
           
  double y24 = y2 - y4;
  double y52 = y5 - y2;
  double y45 = y4 - y5;
  gradop[1][3] = ( x2*(y6-y3-y45) + x3*y24 + x4*(y3-y8-y52)
      + x5*(y8-y6-y24) + x6*y52 + x8*y45 ) / 12.0;
             
  double y31 = y3 - y1;
  double y63 = y6 - y3;
  double y16 = y1 - y6;
  gradop[2][3] = ( x3*(y7-y4-y16) + x4*y31 + x1*(y4-y5-y63)
      + x6*(y5-y7-y31) + x7*y63 + x5*y16 ) / 12.0;
               
  double y42 = y4 - y2;
  double y74 = y7 - y4;
  double y27 = y2 - y7;
  gradop[3][3] = ( x4*(y8-y1-y27) + x1*y42 + x2*(y1-y6-y74)
      + x7*(y6-y8-y42) + x8*y74 + x6*y27 ) / 12.0;
                 
  double y13 = y1 - y3;
  double y81 = y8 - y1;
  double y38 = y3 - y8;
  gradop[4][3] = ( x1*(y5-y2-y38) + x2*y13 + x3*(y2-y7-y81)
      + x8*(y7-y5-y13) + x5*y81 + x7*y38 ) / 12.0;
                   
  double y86 = y8 - y6;
  double y18 = y1 - y8;
  double y61 = y6 - y1;
  gradop[5][3] = ( x8*(y4-y7-y61) + x7*y86 + x6*(y7-y2-y18)
      + x1*(y2-y4-y86) + x4*y18 + x2*y61 ) / 12.0;
                     
  double y57 = y5 - y7;
  double y25 = y2 - y5;
  double y72 = y7 - y2;
  gradop[6][3] = ( x5*(y1-y8-y72) + x8*y57 + x7*(y8-y3-y25)
      + x2*(y3-y1-y57) + x1*y25 + x3*y72 ) / 12.0;
                       
  double y68 = y6 - y8;
  double y36 = y3 - y6;
  double y83 = y8 - y3;
  gradop[7][3] = ( x6*(y2-y5-y83) + x5*y68 + x8*(y5-y4-y36)
      + x3*(y4-y2-y68) + x2*y36 + x4*y83 ) / 12.0;
                         
  double y75 = y7 - y5;
  double y47 = y4 - y7;
  double y54 = y5 - y4;
  gradop[8][3] = ( x7*(y3-y6-y54) + x6*y75 + x5*(y6-y1-y47)
      + x4*(y1-y3-y75) + x3*y47 + x1*y54 ) / 12.0;

  //     calculate element volume and characteristic element aspect ratio
  //     (used in time step and hourglass control) - 



  double volume =  coordinates[0][0] * gradop[1][1]
    + coordinates[1][0] * gradop[2][1]
    + coordinates[2][0] * gradop[3][1]
    + coordinates[3][0] * gradop[4][1]
    + coordinates[4][0] * gradop[5][1]
    + coordinates[5][0] * gradop[6][1]
    + coordinates[6][0] * gradop[7][1]
    + coordinates[7][0] * gradop[8][1];
  double aspect = .5*SQR(volume) /
    ( SQR(gradop[1][1]) + SQR(gradop[2][1])
      + SQR(gradop[3][1]) + SQR(gradop[4][1])
      + SQR(gradop[5][1]) + SQR(gradop[6][1])
      + SQR(gradop[7][1]) + SQR(gradop[8][1])
      + SQR(gradop[1][2]) + SQR(gradop[2][2])
      + SQR(gradop[3][2]) + SQR(gradop[4][2])
      + SQR(gradop[5][2]) + SQR(gradop[6][2])
      + SQR(gradop[7][2]) + SQR(gradop[8][2])
      + SQR(gradop[1][3]) + SQR(gradop[2][3])
      + SQR(gradop[3][3]) + SQR(gradop[4][3])
      + SQR(gradop[5][3]) + SQR(gradop[6][3])
      + SQR(gradop[7][3]) + SQR(gradop[8][3]) );
     
  return (VERDICT_REAL)sqrt(aspect);
  
}

/*!
  oddy of a hex

  General distortion measure based on left Cauchy-Green Tensor
*/
C_FUNC_DEF VERDICT_REAL v_hex_oddy( int /*num_nodes*/, VERDICT_REAL coordinates[][3] )
{
  
  double oddy = 0.0, current_oddy = 0.0;
  VerdictVector xxi, xet, xze;
  
  VerdictVector node_pos[8];
  make_hex_nodes ( coordinates, node_pos );
  
  xxi = calc_hex_efg(1, node_pos);
  xet = calc_hex_efg(2, node_pos);
  xze = calc_hex_efg(3, node_pos);
  
  current_oddy = oddy_comp( xxi, xet, xze);
  if ( current_oddy > oddy ) { oddy = current_oddy; }
  
  xxi.set(coordinates[1][0] - coordinates[0][0],
          coordinates[1][1] - coordinates[0][1],
          coordinates[1][2] - coordinates[0][2] );
  
  xet.set(coordinates[3][0] - coordinates[0][0],
          coordinates[3][1] - coordinates[0][1],
          coordinates[3][2] - coordinates[0][2] );
  
  xze.set(coordinates[4][0] - coordinates[0][0],
          coordinates[4][1] - coordinates[0][1],
          coordinates[4][2] - coordinates[0][2] );
  
  
  current_oddy = oddy_comp( xxi, xet, xze);
  if ( current_oddy > oddy ) { oddy = current_oddy; }
  
  xxi.set(coordinates[2][0] - coordinates[1][0],
          coordinates[2][1] - coordinates[1][1],
          coordinates[2][2] - coordinates[1][2] );
  
  xet.set(coordinates[0][0] - coordinates[1][0],
          coordinates[0][1] - coordinates[1][1],
          coordinates[0][2] - coordinates[1][2] );
  
  xze.set(coordinates[5][0] - coordinates[1][0],
          coordinates[5][1] - coordinates[1][1],
          coordinates[5][2] - coordinates[1][2] );
  
  
  current_oddy = oddy_comp( xxi, xet, xze);
  if ( current_oddy > oddy ) { oddy = current_oddy; }
  
  xxi.set(coordinates[3][0] - coordinates[2][0],
          coordinates[3][1] - coordinates[2][1],
          coordinates[3][2] - coordinates[2][2] );
  
  xet.set(coordinates[1][0] - coordinates[2][0],
          coordinates[1][1] - coordinates[2][1],
          coordinates[1][2] - coordinates[2][2] );
  
  xze.set(coordinates[6][0] - coordinates[2][0],
          coordinates[6][1] - coordinates[2][1],
          coordinates[6][2] - coordinates[2][2] );
  
  
  current_oddy = oddy_comp( xxi, xet, xze);
  if ( current_oddy > oddy ) { oddy = current_oddy; }
  
  xxi.set(coordinates[0][0] - coordinates[3][0],
          coordinates[0][1] - coordinates[3][1],
          coordinates[0][2] - coordinates[3][2] );
  
  xet.set(coordinates[2][0] - coordinates[3][0],
          coordinates[2][1] - coordinates[3][1],
          coordinates[2][2] - coordinates[3][2] );
  
  xze.set(coordinates[7][0] - coordinates[3][0],
          coordinates[7][1] - coordinates[3][1],
          coordinates[7][2] - coordinates[3][2] );
 
 
  current_oddy = oddy_comp( xxi, xet, xze);
  if ( current_oddy > oddy ) { oddy = current_oddy; }
  
  xxi.set(coordinates[7][0] - coordinates[4][0],
          coordinates[7][1] - coordinates[4][1],
          coordinates[7][2] - coordinates[4][2] );
  
  xet.set(coordinates[5][0] - coordinates[4][0],
          coordinates[5][1] - coordinates[4][1],
          coordinates[5][2] - coordinates[4][2] );
  
  xze.set(coordinates[0][0] - coordinates[4][0],
          coordinates[0][1] - coordinates[4][1],
          coordinates[0][2] - coordinates[4][2] );
  
  
  current_oddy = oddy_comp( xxi, xet, xze);
  if ( current_oddy > oddy ) { oddy = current_oddy; }
  
  xxi.set(coordinates[4][0] - coordinates[5][0],
          coordinates[4][1] - coordinates[5][1],
          coordinates[4][2] - coordinates[5][2] );
  
  xet.set(coordinates[6][0] - coordinates[5][0],
          coordinates[6][1] - coordinates[5][1],
          coordinates[6][2] - coordinates[5][2] );
  
  xze.set(coordinates[1][0] - coordinates[5][0],
          coordinates[1][1] - coordinates[5][1],
          coordinates[1][2] - coordinates[5][2] );
  
  
  current_oddy = oddy_comp( xxi, xet, xze);
  if ( current_oddy > oddy ) { oddy = current_oddy; }
  
  xxi.set(coordinates[5][0] - coordinates[6][0],
          coordinates[5][1] - coordinates[6][1],
          coordinates[5][2] - coordinates[6][2] );
  
  xet.set(coordinates[7][0] - coordinates[6][0],
          coordinates[7][1] - coordinates[6][1],
          coordinates[7][2] - coordinates[6][2] );
  
  xze.set(coordinates[2][0] - coordinates[6][0],
          coordinates[2][1] - coordinates[6][1],
          coordinates[2][2] - coordinates[6][2] );
  
  
  current_oddy = oddy_comp( xxi, xet, xze);
  if ( current_oddy > oddy ) { oddy = current_oddy; }
  
  xxi.set(coordinates[6][0] - coordinates[7][0],
          coordinates[6][1] - coordinates[7][1],
          coordinates[6][2] - coordinates[7][2] );
  
  xet.set(coordinates[4][0] - coordinates[7][0],
          coordinates[4][1] - coordinates[7][1],
          coordinates[4][2] - coordinates[7][2] );
  
  xze.set(coordinates[3][0] - coordinates[7][0],
          coordinates[3][1] - coordinates[7][1],
          coordinates[3][2] - coordinates[7][2] );
  
  
  current_oddy = oddy_comp( xxi, xet, xze);
  if ( current_oddy > oddy ) { oddy = current_oddy; }
  
  if( oddy > 0 )
    return (VERDICT_REAL) std::min( oddy, VERDICT_DBL_MAX );
  return (VERDICT_REAL) std::max( oddy, -VERDICT_DBL_MAX );

}

/*!
  condition number of a hex

  Maximum condition number of the Jacobian matrix at 8 corners
*/
C_FUNC_DEF VERDICT_REAL v_hex_condition( int /*num_nodes*/, VERDICT_REAL coordinates[][3] )
{

  VerdictVector node_pos[8];
  make_hex_nodes ( coordinates, node_pos );

  double condition = 0.0, current_condition = 0.0; 
  VerdictVector xxi, xet, xze;

  xxi = calc_hex_efg(1, node_pos );
  xet = calc_hex_efg(2, node_pos );
  xze = calc_hex_efg(3, node_pos );

  current_condition = condition_comp( xxi, xet, xze );
  if ( current_condition > condition ) { condition = current_condition; }

        // J(0,0,0):

  xxi = node_pos[1] - node_pos[0];
  xet = node_pos[3] - node_pos[0];
  xze = node_pos[4] - node_pos[0];

  current_condition = condition_comp( xxi, xet, xze );
  if ( current_condition > condition ) { condition = current_condition; }

        // J(1,0,0):
  
  xxi = node_pos[2] - node_pos[1];
  xet = node_pos[0] - node_pos[1];
  xze = node_pos[5] - node_pos[1];

  current_condition = condition_comp( xxi, xet, xze );
  if ( current_condition > condition ) { condition = current_condition; }

  // J(1,1,0):
  
  xxi = node_pos[3] - node_pos[2];
  xet = node_pos[1] - node_pos[2];
  xze = node_pos[6] - node_pos[2];

  current_condition = condition_comp( xxi, xet, xze );
  if ( current_condition > condition ) { condition = current_condition; }

  // J(0,1,0):
  
  xxi = node_pos[0] - node_pos[3];
  xet = node_pos[2] - node_pos[3];
  xze = node_pos[7] - node_pos[3];

  current_condition = condition_comp( xxi, xet, xze );
  if ( current_condition > condition ) { condition = current_condition; }

  // J(0,0,1):

  xxi = node_pos[7] - node_pos[4];
  xet = node_pos[5] - node_pos[4];
  xze = node_pos[0] - node_pos[4];

  current_condition = condition_comp( xxi, xet, xze );
  if ( current_condition > condition ) { condition = current_condition; }

  // J(1,0,1):

  xxi = node_pos[4] - node_pos[5];
  xet = node_pos[6] - node_pos[5];
  xze = node_pos[1] - node_pos[5];

  current_condition = condition_comp( xxi, xet, xze );
  if ( current_condition > condition ) { condition = current_condition; }

  // J(1,1,1):

  xxi = node_pos[5] - node_pos[6];
  xet = node_pos[7] - node_pos[6];
  xze = node_pos[2] - node_pos[6];

  current_condition = condition_comp( xxi, xet, xze );
  if ( current_condition > condition ) { condition = current_condition; }

  // J(1,1,1):

  xxi = node_pos[6] - node_pos[7];
  xet = node_pos[4] - node_pos[7];
  xze = node_pos[3] - node_pos[7];

  current_condition = condition_comp( xxi, xet, xze );
  if ( current_condition > condition ) { condition = current_condition; }

  condition /= 3.0;

  if( condition > 0 )
    return (VERDICT_REAL) std::min( condition, VERDICT_DBL_MAX );
  return (VERDICT_REAL) std::max( condition, -VERDICT_DBL_MAX );

}

/*!
  jacobian of a hex

  Minimum pointwise volume of local map at 8 corners & center of hex
*/
C_FUNC_DEF VERDICT_REAL v_hex_jacobian( int /*num_nodes*/, VERDICT_REAL coordinates[][3] )
{
  
  VerdictVector node_pos[8];
  make_hex_nodes ( coordinates, node_pos );

  double jacobian = VERDICT_DBL_MAX;
  double current_jacobian = VERDICT_DBL_MAX; 
  VerdictVector xxi, xet, xze;

  xxi = calc_hex_efg(1, node_pos );
  xet = calc_hex_efg(2, node_pos );
  xze = calc_hex_efg(3, node_pos );


  current_jacobian = xxi % (xet * xze) / 64.0;
  if ( current_jacobian < jacobian ) { jacobian = current_jacobian; }

  // J(0,0,0):

  xxi = node_pos[1] - node_pos[0];
  xet = node_pos[3] - node_pos[0];
  xze = node_pos[4] - node_pos[0];

  current_jacobian = xxi % (xet * xze);
  if ( current_jacobian < jacobian ) { jacobian = current_jacobian; }

  // J(1,0,0):
  
  xxi = node_pos[2] - node_pos[1];
  xet = node_pos[0] - node_pos[1];
  xze = node_pos[5] - node_pos[1];

  current_jacobian = xxi % (xet * xze);
  if ( current_jacobian < jacobian ) { jacobian = current_jacobian; }

  // J(1,1,0):
  
  xxi = node_pos[3] - node_pos[2];
  xet = node_pos[1] - node_pos[2];
  xze = node_pos[6] - node_pos[2];

  current_jacobian = xxi % (xet * xze);
  if ( current_jacobian < jacobian ) { jacobian = current_jacobian; }

  // J(0,1,0):
  
  xxi = node_pos[0] - node_pos[3];
  xet = node_pos[2] - node_pos[3];
  xze = node_pos[7] - node_pos[3];

  current_jacobian = xxi % (xet * xze);
  if ( current_jacobian < jacobian ) { jacobian = current_jacobian; }

  // J(0,0,1):

  xxi = node_pos[7] - node_pos[4];
  xet = node_pos[5] - node_pos[4];
  xze = node_pos[0] - node_pos[4];

  current_jacobian = xxi % (xet * xze);
  if ( current_jacobian < jacobian ) { jacobian = current_jacobian; }

  // J(1,0,1):

  xxi = node_pos[4] - node_pos[5];
  xet = node_pos[6] - node_pos[5];
  xze = node_pos[1] - node_pos[5];

  current_jacobian = xxi % (xet * xze);
  if ( current_jacobian < jacobian ) { jacobian = current_jacobian; }

  // J(1,1,1):

  xxi = node_pos[5] - node_pos[6];
  xet = node_pos[7] - node_pos[6];
  xze = node_pos[2] - node_pos[6];

  current_jacobian = xxi % (xet * xze);
  if ( current_jacobian < jacobian ) { jacobian = current_jacobian; }

  // J(0,1,1):

  xxi = node_pos[6] - node_pos[7];
  xet = node_pos[4] - node_pos[7];
  xze = node_pos[3] - node_pos[7];

  current_jacobian = xxi % (xet * xze);
  if ( current_jacobian < jacobian ) { jacobian = current_jacobian; }

  if( jacobian > 0 )
    return (VERDICT_REAL) std::min( jacobian, VERDICT_DBL_MAX );
  return (VERDICT_REAL) std::max( jacobian, -VERDICT_DBL_MAX );
}

/*!
  scaled jacobian of a hex

  Minimum Jacobian divided by the lengths of the 3 edge vectors
*/
C_FUNC_DEF VERDICT_REAL v_hex_scaled_jacobian( int /*num_nodes*/, VERDICT_REAL coordinates[][3] )
{

  double jacobi, min_norm_jac = VERDICT_DBL_MAX;  
  double min_jacobi = VERDICT_DBL_MAX;
  double temp_norm_jac, lengths;
  double len1_sq, len2_sq, len3_sq; 
  VerdictVector xxi, xet, xze;

  VerdictVector node_pos[8];
  make_hex_nodes ( coordinates, node_pos );

  xxi = calc_hex_efg(1, node_pos );
  xet = calc_hex_efg(2, node_pos );
  xze = calc_hex_efg(3, node_pos );

  jacobi = xxi % ( xet * xze );
  if( jacobi < min_jacobi) { min_jacobi = jacobi; }

  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return (VERDICT_REAL) VERDICT_DBL_MAX;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  temp_norm_jac = jacobi / lengths;
  
  if( temp_norm_jac < min_norm_jac)
    min_norm_jac = temp_norm_jac;  
  else  
    temp_norm_jac = jacobi;

  // J(0,0,0):

  xxi = node_pos[1] - node_pos[0];
  xet = node_pos[3] - node_pos[0];
  xze = node_pos[4] - node_pos[0];

  jacobi = xxi % ( xet * xze );
  if( jacobi < min_jacobi ) { min_jacobi = jacobi; }

  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return (VERDICT_REAL) VERDICT_DBL_MAX;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  temp_norm_jac = jacobi / lengths;
  if( temp_norm_jac < min_norm_jac)
    min_norm_jac = temp_norm_jac;  
  else  
    temp_norm_jac = jacobi;


  // J(1,0,0):
  
  xxi = node_pos[2] - node_pos[1];
  xet = node_pos[0] - node_pos[1];
  xze = node_pos[5] - node_pos[1];

  jacobi = xxi % ( xet * xze );
  if( jacobi < min_jacobi ) { min_jacobi = jacobi; }
  
  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return (VERDICT_REAL) VERDICT_DBL_MAX;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  temp_norm_jac = jacobi / lengths;
  if( temp_norm_jac < min_norm_jac)
    min_norm_jac = temp_norm_jac;  
  else  
    temp_norm_jac = jacobi;

  // J(1,1,0):
  
  xxi = node_pos[3] - node_pos[2];
  xet = node_pos[1] - node_pos[2];
  xze = node_pos[6] - node_pos[2];

  jacobi = xxi % ( xet * xze );
  if( jacobi < min_jacobi ) { min_jacobi = jacobi; }
  
  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return (VERDICT_REAL) VERDICT_DBL_MAX;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  temp_norm_jac = jacobi / lengths;
  if( temp_norm_jac < min_norm_jac)
    min_norm_jac = temp_norm_jac;  
  else  
    temp_norm_jac = jacobi;

  // J(0,1,0):
  
  xxi = node_pos[0] - node_pos[3];
  xet = node_pos[2] - node_pos[3];
  xze = node_pos[7] - node_pos[3];

  jacobi = xxi % ( xet * xze );
  if( jacobi < min_jacobi ) { min_jacobi = jacobi; }
  
  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return (VERDICT_REAL) VERDICT_DBL_MAX;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  temp_norm_jac = jacobi / lengths;
  if( temp_norm_jac < min_norm_jac)
    min_norm_jac = temp_norm_jac;  
  else  
    temp_norm_jac = jacobi;

  // J(0,0,1):

  xxi = node_pos[7] - node_pos[4];
  xet = node_pos[5] - node_pos[4];
  xze = node_pos[0] - node_pos[4];

  jacobi = xxi % ( xet * xze );
  if( jacobi < min_jacobi ) { min_jacobi = jacobi; }
  
  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return (VERDICT_REAL) VERDICT_DBL_MAX;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  temp_norm_jac = jacobi / lengths;
  if( temp_norm_jac < min_norm_jac)
    min_norm_jac = temp_norm_jac;  
  else 
    temp_norm_jac = jacobi;

  // J(1,0,1):

  xxi = node_pos[4] - node_pos[5];
  xet = node_pos[6] - node_pos[5];
  xze = node_pos[1] - node_pos[5];

  jacobi = xxi % ( xet * xze );
  if( jacobi < min_jacobi ) { min_jacobi = jacobi; }
  
  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return (VERDICT_REAL) VERDICT_DBL_MAX;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  temp_norm_jac = jacobi / lengths;
  if( temp_norm_jac < min_norm_jac)
    min_norm_jac = temp_norm_jac;  
  else
    temp_norm_jac = jacobi;

  // J(1,1,1):

  xxi = node_pos[5] - node_pos[6];
  xet = node_pos[7] - node_pos[6];
  xze = node_pos[2] - node_pos[6];

  jacobi = xxi % ( xet * xze );
  if( jacobi < min_jacobi ) { min_jacobi = jacobi; }
  
  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return (VERDICT_REAL) VERDICT_DBL_MAX;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  temp_norm_jac = jacobi / lengths;
  if( temp_norm_jac < min_norm_jac)
    min_norm_jac = temp_norm_jac;  
  else
    temp_norm_jac = jacobi;

  // J(0,1,1):

  xxi = node_pos[6] - node_pos[7];
  xet = node_pos[4] - node_pos[7];
  xze = node_pos[3] - node_pos[7];

  jacobi = xxi % ( xet * xze );
  if( jacobi < min_jacobi ) { min_jacobi = jacobi; }
  
  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return (VERDICT_REAL) VERDICT_DBL_MAX;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  temp_norm_jac = jacobi / lengths;
  if( temp_norm_jac < min_norm_jac)
    min_norm_jac = temp_norm_jac;  
  else 
    temp_norm_jac = jacobi;

  if( min_norm_jac> 0 )
    return (VERDICT_REAL) std::min( min_norm_jac, VERDICT_DBL_MAX );
  return (VERDICT_REAL) std::max( min_norm_jac, -VERDICT_DBL_MAX );
}

/*!
  shear of a hex
  
  3/Condition number of Jacobian Skew matrix
*/
C_FUNC_DEF VERDICT_REAL v_hex_shear( int /*num_nodes*/, VERDICT_REAL coordinates[][3] )
{

  double shear;
  double min_shear = 1.0; 
  VerdictVector xxi, xet, xze;
  double det, len1_sq, len2_sq, len3_sq, lengths;


  VerdictVector node_pos[8];
  make_hex_nodes ( coordinates, node_pos );

  // J(0,0,0):

  xxi = node_pos[1] - node_pos[0];
  xet = node_pos[3] - node_pos[0];
  xze = node_pos[4] - node_pos[0];

  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return 0;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  det = xxi % (xet * xze);
  if ( det < VERDICT_DBL_MIN ) { return 0; }

  shear = det / lengths;
  min_shear = std::min( shear, min_shear );


  // J(1,0,0):
  
  xxi = node_pos[2] - node_pos[1];
  xet = node_pos[0] - node_pos[1];
  xze = node_pos[5] - node_pos[1];

  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return 0;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  det = xxi % (xet * xze);
  if ( det < VERDICT_DBL_MIN ) { return 0; }

  shear = det / lengths; 
  min_shear = std::min( shear, min_shear );


  // J(1,1,0):
  
  xxi = node_pos[3] - node_pos[2];
  xet = node_pos[1] - node_pos[2];
  xze = node_pos[6] - node_pos[2];

  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return 0;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  det = xxi % (xet * xze);
  if ( det < VERDICT_DBL_MIN ) { return 0; }

  shear = det / lengths; 
  min_shear = std::min( shear, min_shear );


  // J(0,1,0):
  
  xxi = node_pos[0] - node_pos[3];
  xet = node_pos[2] - node_pos[3];
  xze = node_pos[7] - node_pos[3];

  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return 0;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  det = xxi % (xet * xze);
  if ( det < VERDICT_DBL_MIN ) { return 0; }

  shear = det / lengths; 
  min_shear = std::min( shear, min_shear );


  // J(0,0,1):

  xxi = node_pos[7] - node_pos[4];
  xet = node_pos[5] - node_pos[4];
  xze = node_pos[0] - node_pos[4];

  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return 0;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  det = xxi % (xet * xze);
  if ( det < VERDICT_DBL_MIN ) { return 0; }

  shear = det / lengths; 
  min_shear = std::min( shear, min_shear );


  // J(1,0,1):

  xxi = node_pos[4] - node_pos[5];
  xet = node_pos[6] - node_pos[5];
  xze = node_pos[1] - node_pos[5];

  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return 0;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  det = xxi % (xet * xze);
  if ( det < VERDICT_DBL_MIN ) { return 0; }

  shear = det / lengths; 
  min_shear = std::min( shear, min_shear );


  // J(1,1,1):

  xxi = node_pos[5] - node_pos[6];
  xet = node_pos[7] - node_pos[6];
  xze = node_pos[2] - node_pos[6];

  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return 0;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  det = xxi % (xet * xze);
  if ( det < VERDICT_DBL_MIN ) { return 0; }

  shear = det / lengths; 
  min_shear = std::min( shear, min_shear );


  // J(0,1,1):

  xxi = node_pos[6] - node_pos[7];
  xet = node_pos[4] - node_pos[7];
  xze = node_pos[3] - node_pos[7];

  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return 0;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  det = xxi % (xet * xze);
  if ( det < VERDICT_DBL_MIN ) { return 0; }

  shear = det / lengths; 
  min_shear = std::min( shear, min_shear );

  if( min_shear <= VERDICT_DBL_MIN )
    min_shear = 0;
  
  if( min_shear > 0 )
    return (VERDICT_REAL) std::min( min_shear, VERDICT_DBL_MAX );
  return (VERDICT_REAL) std::max( min_shear, -VERDICT_DBL_MAX );
}

/*!
  shape of a hex

  3/Condition number of weighted Jacobian matrix
*/
C_FUNC_DEF VERDICT_REAL v_hex_shape( int /*num_nodes*/, VERDICT_REAL coordinates[][3] )
{


  double det, shape;
  double min_shape = 1.0; 
  static const double two_thirds = 2.0/3.0;

  VerdictVector xxi, xet, xze;

  VerdictVector node_pos[8];
  make_hex_nodes ( coordinates, node_pos );


  // J(0,0,0):

  xxi = node_pos[1] - node_pos[0];
  xet = node_pos[3] - node_pos[0];
  xze = node_pos[4] - node_pos[0];

  det = xxi % (xet * xze);
  if ( det > VERDICT_DBL_MIN ) 
    shape = 3 * pow( det, two_thirds) / (xxi%xxi + xet%xet + xze%xze);
  else
    return 0;
  
  if( shape < min_shape ) { min_shape = shape; }
  

  // J(1,0,0):
  
  xxi = node_pos[2] - node_pos[1];
  xet = node_pos[0] - node_pos[1];
  xze = node_pos[5] - node_pos[1];

  det = xxi % (xet * xze);
  if ( det > VERDICT_DBL_MIN ) 
    shape = 3 * pow( det, two_thirds) / (xxi%xxi + xet%xet + xze%xze);
  else
    return 0;
  
  if( shape < min_shape ) { min_shape = shape; }
  

  // J(1,1,0):
  
  xxi = node_pos[3] - node_pos[2];
  xet = node_pos[1] - node_pos[2];
  xze = node_pos[6] - node_pos[2];

  det = xxi % (xet * xze);
  if ( det > VERDICT_DBL_MIN ) 
    shape = 3 * pow( det, two_thirds) / (xxi%xxi + xet%xet + xze%xze);
  else
    return 0;
  
  if( shape < min_shape ) { min_shape = shape; }
  

  // J(0,1,0):
  
  xxi = node_pos[0] - node_pos[3];
  xet = node_pos[2] - node_pos[3];
  xze = node_pos[7] - node_pos[3];

  det = xxi % (xet * xze);
  if ( det > VERDICT_DBL_MIN ) 
    shape = 3 * pow( det, two_thirds) / (xxi%xxi + xet%xet + xze%xze);
  else
    return 0;
  
  if( shape < min_shape ) { min_shape = shape; }
  

  // J(0,0,1):

  xxi = node_pos[7] - node_pos[4];
  xet = node_pos[5] - node_pos[4];
  xze = node_pos[0] - node_pos[4];

  det = xxi % (xet * xze);
  if ( det > VERDICT_DBL_MIN ) 
    shape = 3 * pow( det, two_thirds) / (xxi%xxi + xet%xet + xze%xze);
  else
    return 0;
  
  if( shape < min_shape ) { min_shape = shape; }
  

  // J(1,0,1):

  xxi = node_pos[4] - node_pos[5];
  xet = node_pos[6] - node_pos[5];
  xze = node_pos[1] - node_pos[5];

  det = xxi % (xet * xze);
  if ( det > VERDICT_DBL_MIN ) 
    shape = 3 * pow( det, two_thirds) / (xxi%xxi + xet%xet + xze%xze);
  else
    return 0;
  
  if( shape < min_shape ) { min_shape = shape; }
  

  // J(1,1,1):

  xxi = node_pos[5] - node_pos[6];
  xet = node_pos[7] - node_pos[6];
  xze = node_pos[2] - node_pos[6];

  det = xxi % (xet * xze);
  if ( det > VERDICT_DBL_MIN ) 
    shape = 3 * pow( det, two_thirds) / (xxi%xxi + xet%xet + xze%xze);
  else
    return 0;
  
  if( shape < min_shape ) { min_shape = shape; }
  

  // J(1,1,1):

  xxi = node_pos[6] - node_pos[7];
  xet = node_pos[4] - node_pos[7];
  xze = node_pos[3] - node_pos[7];

  det = xxi % (xet * xze);
  if ( det > VERDICT_DBL_MIN ) 
    shape = 3 * pow( det, two_thirds) / (xxi%xxi + xet%xet + xze%xze);
  else
    return 0;
  
  if( shape < min_shape ) { min_shape = shape; }

  if( min_shape <= VERDICT_DBL_MIN )
    min_shape = 0;

  if( min_shape > 0 )
    return (VERDICT_REAL) std::min( min_shape, VERDICT_DBL_MAX );
  return (VERDICT_REAL) std::max( min_shape, -VERDICT_DBL_MAX );
}

/*!
  relative size of a hex

  Min( J, 1/J ), where J is determinant of weighted Jacobian matrix
*/
C_FUNC_DEF VERDICT_REAL v_hex_relative_size_squared( int /*num_nodes*/, VERDICT_REAL coordinates[][3] )
{
  double size = 0;
  double tau; 

  VerdictVector xxi, xet, xze;
  double det, det_sum = 0;

  v_hex_get_weight( xxi, xet, xze );
 
  //This is the average relative size 
  double detw = xxi % (xet * xze);

  if ( detw < VERDICT_DBL_MIN ) 
    return 0; 

  VerdictVector node_pos[8];
  make_hex_nodes ( coordinates, node_pos );

  // J(0,0,0):

  xxi = node_pos[1] - node_pos[0];
  xet = node_pos[3] - node_pos[0];
  xze = node_pos[4] - node_pos[0];
  
  det = xxi % (xet * xze);
  det_sum += det;  


  // J(1,0,0):
  
  xxi = node_pos[2] - node_pos[1];
  xet = node_pos[0] - node_pos[1];
  xze = node_pos[5] - node_pos[1];

  det = xxi % (xet * xze);
  det_sum += det;  


  // J(0,1,0):
  
  xxi = node_pos[3] - node_pos[2];
  xet = node_pos[1] - node_pos[2];
  xze = node_pos[6] - node_pos[2];

  det = xxi % (xet * xze);
  det_sum += det;  

  // J(1,1,0):
  
  xxi = node_pos[0] - node_pos[3];
  xet = node_pos[2] - node_pos[3];
  xze = node_pos[7] - node_pos[3];

  det = xxi % (xet * xze);
  det_sum += det;  


  // J(0,1,0):

  xxi = node_pos[7] - node_pos[4];
  xet = node_pos[5] - node_pos[4];
  xze = node_pos[0] - node_pos[4];

  det = xxi % (xet * xze);
  det_sum += det;  


  // J(1,0,1):

  xxi = node_pos[4] - node_pos[5];
  xet = node_pos[6] - node_pos[5];
  xze = node_pos[1] - node_pos[5];

  det = xxi % (xet * xze);
  det_sum += det;  


  // J(1,1,1):

  xxi = node_pos[5] - node_pos[6];
  xet = node_pos[7] - node_pos[6];
  xze = node_pos[2] - node_pos[6];

  det = xxi % (xet * xze);
  det_sum += det;  


  // J(1,1,1):

  xxi = node_pos[6] - node_pos[7];
  xet = node_pos[4] - node_pos[7];
  xze = node_pos[3] - node_pos[7];

  det = xxi % (xet * xze);
  det_sum += det;  


  if( det_sum > VERDICT_DBL_MIN )
  {
    tau = det_sum / ( 8*detw);

    tau = std::min( tau, 1.0/tau); 

    size = tau*tau; 
  }

  if( size > 0 )
    return (VERDICT_REAL) std::min( size, VERDICT_DBL_MAX );
  return (VERDICT_REAL) std::max( size, -VERDICT_DBL_MAX );
}

/*!
  shape and size of a hex

  Product of Shape and Relative Size
*/
C_FUNC_DEF VERDICT_REAL v_hex_shape_and_size( int num_nodes, VERDICT_REAL coordinates[][3] )
{
  double size = v_hex_relative_size_squared( num_nodes, coordinates );
  double shape = v_hex_shape( num_nodes, coordinates );

  double shape_size = size * shape;

  if( shape_size > 0 )
    return (VERDICT_REAL) std::min( shape_size, VERDICT_DBL_MAX );
  return (VERDICT_REAL) std::max( shape_size, -VERDICT_DBL_MAX );

}



/*!
  shear and size of a hex

  Product of Shear and Relative Size
*/
C_FUNC_DEF VERDICT_REAL v_hex_shear_and_size( int num_nodes, VERDICT_REAL coordinates[][3] )
{
  double size = v_hex_relative_size_squared( num_nodes, coordinates );
  double shear = v_hex_shear( num_nodes, coordinates );

  double shear_size = shear * size; 

  if( shear_size > 0 )
    return (VERDICT_REAL) std::min( shear_size, VERDICT_DBL_MAX );
  return (VERDICT_REAL) std::max( shear_size, -VERDICT_DBL_MAX );

}

/*!
  distortion of a hex
*/
VERDICT_REAL v_hex_distortion( int num_nodes, VERDICT_REAL coordinates[][3] )
{

   //use 2x2 gauss points for linear hex and 3x3 for 2nd order hex
   int number_of_gauss_points=0;
   if (num_nodes ==8)
      //2x2 quadrature rule
      number_of_gauss_points = 2;
   else if (num_nodes ==20)
      //3x3 quadrature rule
      number_of_gauss_points = 3;

   int number_dimension = 3;
   int total_number_of_gauss_points = number_of_gauss_points
      *number_of_gauss_points*number_of_gauss_points;
   double distortion = VERDICT_DBL_MAX;

   // maxTotalNumberGaussPoints =27, maxNumberNodes = 20
   // they are defined in GaussIntegration.h
   // This is used to make these arrays static.
   // I tried dynamically allocated arrays but the new and delete
   // was very expensive

   double shape_function[maxTotalNumberGaussPoints][maxNumberNodes];
   double dndy1[maxTotalNumberGaussPoints][maxNumberNodes];
   double dndy2[maxTotalNumberGaussPoints][maxNumberNodes];
   double dndy3[maxTotalNumberGaussPoints][maxNumberNodes];
   double weight[maxTotalNumberGaussPoints];


   //create an object of GaussIntegration
   GaussIntegration::initialize(number_of_gauss_points,num_nodes,number_dimension );
   GaussIntegration::calculate_shape_function_3d_hex();
   GaussIntegration::get_shape_func(shape_function[0], dndy1[0], dndy2[0], dndy3[0],weight);


   VerdictVector xxi, xet, xze, xin;

   double jacobian, minimum_jacobian;
   double element_volume =0.0;
   minimum_jacobian = VERDICT_DBL_MAX;
   // calculate element volume
   int ife, ja;
   for (ife=0;ife<total_number_of_gauss_points; ife++)
   {

      xxi.set(0.0,0.0,0.0);
      xet.set(0.0,0.0,0.0);
      xze.set(0.0,0.0,0.0);

      for (ja=0;ja<num_nodes;ja++)
      {
   xin.set(coordinates[ja][0], coordinates[ja][1], coordinates[ja][2]);
         xxi += dndy1[ife][ja]*xin;
         xet += dndy2[ife][ja]*xin;
         xze += dndy3[ife][ja]*xin;
      }

      jacobian = xxi % (xet * xze);
      if (minimum_jacobian > jacobian)
         minimum_jacobian = jacobian;

      element_volume += weight[ife]*jacobian;
      }

   // loop through all nodes
   double dndy1_at_node[maxNumberNodes][maxNumberNodes];
   double dndy2_at_node[maxNumberNodes][maxNumberNodes];
   double dndy3_at_node[maxNumberNodes][maxNumberNodes];

   GaussIntegration::calculate_derivative_at_nodes_3d( dndy1_at_node, dndy2_at_node, dndy3_at_node);
   int node_id;
   for (node_id=0;node_id<num_nodes; node_id++)
   {

      xxi.set(0.0,0.0,0.0);
      xet.set(0.0,0.0,0.0);
      xze.set(0.0,0.0,0.0);

      for (ja=0;ja<num_nodes;ja++)
      {
   xin.set(coordinates[ja][0], coordinates[ja][1], coordinates[ja][2]);
         xxi += dndy1_at_node[node_id][ja]*xin;
         xet += dndy2_at_node[node_id][ja]*xin;
         xze += dndy3_at_node[node_id][ja]*xin;
      }

      jacobian = xxi % (xet * xze);
      if (minimum_jacobian > jacobian)
         minimum_jacobian = jacobian;

      }
   distortion = minimum_jacobian/element_volume*8.;
   return (VERDICT_REAL)distortion;
}





/*
C_FUNC_DEF VERDICT_REAL hex_jac_normjac_oddy_cond( int choices[], 
                      VERDICT_REAL coordinates[][3],
                      VERDICT_REAL answers[4]  )
{

  //Define variables
  int i;
  
  VERDICT_REAL xxi[3], xet[3], xze[3];
  VERDICT_REAL norm_jacobian = 0.0, current_norm_jac = 0.0;
        VERDICT_REAL jacobian = 0.0, current_jacobian = 0.0;
  VERDICT_REAL oddy = 0.0, current_oddy = 0.0;  
  VERDICT_REAL condition = 0.0, current_condition = 0.0;


        for( i=0; i<3; i++)
          xxi[i] = calc_hex_efg( 2, i, coordinates );
        for( i=0; i<3; i++)
          xet[i] = calc_hex_efg( 3, i, coordinates );
        for( i=0; i<3; i++)
          xze[i] = calc_hex_efg( 6, i, coordinates );

  // norm jacobian and jacobian
  if( choices[0] || choices[1] )
  {
    current_jacobian = determinant( xxi, xet, xze  );
    current_norm_jac = normalize_jacobian( current_jacobian,
            xxi, xet, xze );
    
    if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 

    current_jacobian /= 64.0;
    if (current_jacobian < jacobian) { jacobian = current_jacobian; }
  }

  // oddy 
  if( choices[2] )
  {
    current_oddy = oddy_comp( xxi, xet, xze );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
  }

  // condition
  if( choices[3] )
  {
    current_condition = condition_comp( xxi, xet, xze );
    if ( current_condition > condition ) { condition = current_condition; }
  }


  for( i=0; i<3; i++)
  {
    xxi[i] = coordinates[1][i] - coordinates[0][i];
    xet[i] = coordinates[3][i] - coordinates[0][i];
    xze[i] = coordinates[4][i] - coordinates[0][i];
  } 

  // norm jacobian and jacobian
  if( choices[0] || choices[1] )
  {
    current_jacobian = determinant( xxi, xet, xze  );
    current_norm_jac = normalize_jacobian( current_jacobian,
            xxi, xet, xze );
    
    if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 
    if (current_jacobian < jacobian) { jacobian = current_jacobian; }
  }

  // oddy 
  if( choices[2] )
  {
    current_oddy = oddy_comp( xxi, xet, xze );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
  }

  // condition
  if( choices[3] )
  {
    current_condition = condition_comp( xxi, xet, xze );
    if ( current_condition > condition ) { condition = current_condition; }
  }

  
  for( i=0; i<3; i++)
  {
          xxi[i] = coordinates[1][i] - coordinates[0][i];
          xet[i] = coordinates[2][i] - coordinates[1][i];
          xze[i] = coordinates[5][i] - coordinates[1][i];
  }

  // norm jacobian and jacobian
  if( choices[0] || choices[1] )
  {
    current_jacobian = determinant( xxi,  xet, xze );
    current_norm_jac = normalize_jacobian( current_jacobian,
            xxi, xet, xze );
    
    if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 
    if (current_jacobian < jacobian) { jacobian = current_jacobian; }
  }

  // oddy 
  if( choices[2] )
  {
    current_oddy = oddy_comp( xxi, xet, xze );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
  }

  // condition
  if( choices[3] )
  {
    current_condition = condition_comp( xxi, xet, xze );
    if ( current_condition > condition ) { condition = current_condition; }
  }


  for( i=0; i<3; i++)
  {
          xxi[i] = coordinates[2][i] - coordinates[3][i];
          xet[i] = coordinates[3][i] - coordinates[0][i];
          xze[i] = coordinates[7][i] - coordinates[3][i];
  }

  // norm jacobian and jacobian
  if( choices[0] || choices[1] )
  {
    current_jacobian = determinant( xxi, xet, xze );
    current_norm_jac = normalize_jacobian( current_jacobian,
            xxi, xet, xze );
    
    if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 
    if (current_jacobian < jacobian) { jacobian = current_jacobian; }
  }

  // oddy 
  if( choices[2] )
  {
    current_oddy = oddy_comp( xxi, xet, xze );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
  }

  // condition
  if( choices[3] )
  {
    current_condition = condition_comp( xxi, xet, xze );
    if ( current_condition > condition ) { condition = current_condition; }
  }

  
  for( i=0; i<3; i++)
  {
          xxi[i] = coordinates[5][i] - coordinates[4][i];
          xet[i] = coordinates[7][i] - coordinates[4][i];
          xze[i] = coordinates[4][i] - coordinates[0][i];
  }  


  // norm jacobian and jacobian
  if( choices[0] || choices[1] )
  {
    current_jacobian = determinant( xxi, xet, xze );
    current_norm_jac = normalize_jacobian( current_jacobian,
            xxi, xet, xze );
    
    if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 
    if (current_jacobian < jacobian) { jacobian = current_jacobian; }
  }

  // oddy 
  if( choices[2] )
  {
    current_oddy = oddy_comp( xxi, xet, xze );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
  }

  // condition
  if( choices[3] )
  {
    current_condition = condition_comp( xxi, xet, xze );
    if ( current_condition > condition ) { condition = current_condition; }
  }

  
  for( i=0; i<3; i++)
  {
          xxi[i] = coordinates[2][i] - coordinates[3][i];
          xet[i] = coordinates[2][i] - coordinates[1][i];
          xze[i] = coordinates[6][i] - coordinates[2][i];
  }  

  // norm jacobian and jacobian
  if( choices[0] || choices[1] )
  {
    current_jacobian = determinant( xxi, xet, xze );
    current_norm_jac = normalize_jacobian( current_jacobian,
            xxi, xet, xze );
    
    if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 
    if (current_jacobian < jacobian) { jacobian = current_jacobian; }
  }

  // oddy 
  if( choices[2] )
  {
    current_oddy = oddy_comp( xxi, xet, xze );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
  }

  // condition
  if( choices[3] )
  {
    current_condition = condition_comp( xxi, xet, xze );
    if ( current_condition > condition ) { condition = current_condition; }
  }


  for( i=0; i<3; i++)
  {
          xxi[i] = coordinates[5][i] - coordinates[4][i];
          xet[i] = coordinates[6][i] - coordinates[5][i];
          xze[i] = coordinates[5][i] - coordinates[1][i];
  }


  // norm jacobian and jacobian
  if( choices[0] || choices[1] )
  {
    current_jacobian = determinant( xxi, xet, xze );
    current_norm_jac = normalize_jacobian( current_jacobian,
            xxi, xet, xze );
    
    if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 
    if (current_jacobian < jacobian) { jacobian = current_jacobian; }
  }

  // oddy 
  if( choices[2] )
  {
    current_oddy = oddy_comp( xxi, xet, xze );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
  }

  // condition
  if( choices[3] )
  {
    current_condition = condition_comp( xxi, xet, xze );
    if ( current_condition > condition ) { condition = current_condition; }
  }

    
  for( i=0; i<3; i++)
  {
          xxi[i] = coordinates[6][i] - coordinates[7][i];
          xet[i] = coordinates[7][i] - coordinates[4][i];
          xze[i] = coordinates[7][i] - coordinates[3][i];
  }


  // norm jacobian and jacobian
  if( choices[0] || choices[1] )
  {
    current_jacobian = determinant( xxi, xet, xze );
    current_norm_jac = normalize_jacobian( current_jacobian,
            xxi, xet, xze );
    
    if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 
    if (current_jacobian < jacobian) { jacobian = current_jacobian; }
  }

  // oddy 
  if( choices[2] )
  {
    current_oddy = oddy_comp( xxi, xet, xze );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
  }

  // condition
  if( choices[3] )
  {
    current_condition = condition_comp( xxi, xet, xze );
    if ( current_condition > condition ) { condition = current_condition; }
  }

  
  for( i=0; i<3; i++)
  {
          xxi[i] = coordinates[6][i] - coordinates[7][i];
          xet[i] = coordinates[6][i] - coordinates[5][i];
          xze[i] = coordinates[6][i] - coordinates[2][i];
  }

  // norm jacobian and jacobian
  if( choices[0] || choices[1] )
  {
    current_jacobian = determinant( xxi, xet, xze );
    current_norm_jac = normalize_jacobian( current_jacobian,
            xxi, xet, xze );
    
    if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 
    if (current_jacobian < jacobian) { jacobian = current_jacobian; }
  }

  // oddy 
  if( choices[2] )
  {
    current_oddy = oddy_comp( xxi, xet, xze );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
  }

  // condition
  if( choices[3] )
  {
    current_condition = condition_comp( xxi, xet, xze );
    if ( current_condition > condition ) { condition = current_condition; }

    condition /= 3.0;
  }

  
  answers[0] = jacobian;
  answers[1] = norm_jacobian;
  answers[2] = oddy;
  answers[3] = condition;

  return 1.0;

}
*/

/*!
  multiple quality metrics of a hex
*/
C_FUNC_DEF void v_hex_quality( int num_nodes, VERDICT_REAL coordinates[][3], 
  unsigned int metrics_request_flag, HexMetricVals *metric_vals )
{ 
  memset( metric_vals, 0, sizeof(HexMetricVals) );

  // aspect, skew, taper, and volume
  if(metrics_request_flag & (V_HEX_ASPECT | V_HEX_SKEW | V_HEX_TAPER ) ) 
  {
    VerdictVector node_pos[8];
    make_hex_nodes ( coordinates, node_pos );
   
    VerdictVector efg1, efg2, efg3;
    efg1 = calc_hex_efg( 1, node_pos);
    efg2 = calc_hex_efg( 2, node_pos);
    efg3 = calc_hex_efg( 3, node_pos);

    if(metrics_request_flag & V_HEX_ASPECT)
    {
      double aspect_12, aspect_13, aspect_23;

      double mag_efg1 = efg1.length();
      double mag_efg2 = efg2.length();
      double mag_efg3 = efg3.length();
      
      aspect_12 = safe_ratio( std::max( mag_efg1, mag_efg2 ) , std::min( mag_efg1, mag_efg2 ) );
      aspect_13 = safe_ratio( std::max( mag_efg1, mag_efg3 ) , std::min( mag_efg1, mag_efg3 ) );
      aspect_23 = safe_ratio( std::max( mag_efg2, mag_efg3 ) , std::min( mag_efg2, mag_efg3 ) );

      metric_vals->aspect = (VERDICT_REAL)std::max( aspect_12, std::max( aspect_13, aspect_23 ) );
    }
    
    if(metrics_request_flag & V_HEX_SKEW)
    {
      
      VerdictVector vec1 = efg1;
      VerdictVector vec2 = efg2;
      VerdictVector vec3 = efg3;

      if( vec1.normalize() <= VERDICT_DBL_MIN ||
          vec2.normalize() <= VERDICT_DBL_MIN ||
          vec3.normalize() <= VERDICT_DBL_MIN  )
      {
        metric_vals->skew = (VERDICT_REAL)VERDICT_DBL_MAX;
      }
      else
      {
        double skewx = fabs(vec1 % vec2);
        double skewy = fabs(vec1 % vec3);
        double skewz = fabs(vec2 % vec3);

        metric_vals->skew = (VERDICT_REAL)(std::max( skewx, std::max( skewy, skewz ) ));
      }
    }
  
    if(metrics_request_flag & V_HEX_TAPER)
    {
      VerdictVector efg12 = calc_hex_efg( 12, node_pos);
      VerdictVector efg13 = calc_hex_efg( 13, node_pos);
      VerdictVector efg23 = calc_hex_efg( 23, node_pos);

      double taperx = fabs( safe_ratio( efg12.length(), std::min( efg1.length(), efg2.length())));
      double tapery = fabs( safe_ratio( efg13.length(), std::min( efg1.length(), efg3.length())));
      double taperz = fabs( safe_ratio( efg23.length(), std::min( efg2.length(), efg3.length())));
      
      metric_vals->taper = (VERDICT_REAL)std::max(taperx, std::max(tapery, taperz));  
    }
  }
  
  if(metrics_request_flag & V_HEX_VOLUME)
  {
    metric_vals->volume = v_hex_volume(8, coordinates ); 
  }

  if(metrics_request_flag & ( V_HEX_JACOBIAN | V_HEX_SCALED_JACOBIAN | 
        V_HEX_CONDITION |V_HEX_SHEAR | V_HEX_SHAPE | V_HEX_RELATIVE_SIZE_SQUARED 
        | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE |
        V_HEX_STRETCH ))
  {

    static const double two_thirds = 2.0/3.0;
    VerdictVector edges[12];
    // the length squares
    double length_squared[12];
    // make vectors from coordinates
    make_hex_edges(coordinates, edges);

    // calculate the length squares if we need to
    if(metrics_request_flag & ( V_HEX_JACOBIAN | V_HEX_SHEAR | V_HEX_SCALED_JACOBIAN | V_HEX_SHAPE | 
          V_HEX_SHAPE_AND_SIZE | V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHEAR_AND_SIZE | V_HEX_STRETCH))
      make_edge_length_squares(edges, length_squared);

    double jacobian = VERDICT_DBL_MAX, scaled_jacobian = VERDICT_DBL_MAX,
      condition = 0.0, shear=1.0, shape=1.0, oddy = 0.0;
    double current_jacobian, current_scaled_jacobian, current_condition, current_shape,
      detw=0, det_sum=0, current_oddy;
    VerdictBoolean rel_size_error = VERDICT_FALSE;
    
    VerdictVector xxi, xet, xze;

    // get weights if we need based on average size of a hex
    if(metrics_request_flag & (V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ))
    {
      v_hex_get_weight(xxi, xet, xze);
      detw = xxi % (xet * xze);
      if(detw < VERDICT_DBL_MIN)
        rel_size_error = VERDICT_TRUE;
    }

    xxi = edges[0] - edges[2] + edges[4] - edges[6];
    xet = edges[1] - edges[3] + edges[5] - edges[7];
    xze = edges[8] + edges[9] + edges[10] + edges[11];

    current_jacobian = xxi % (xet * xze) / 64.0;
    if ( current_jacobian < jacobian ) 
      jacobian = current_jacobian;


    if(metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ))
    {
      current_jacobian *= 64.0;
      current_scaled_jacobian = current_jacobian / 
        sqrt(xxi.length_squared() * xet.length_squared() * xze.length_squared());
      if(current_scaled_jacobian < scaled_jacobian)
        shear = scaled_jacobian = current_scaled_jacobian;
    }

    if(metrics_request_flag & V_HEX_CONDITION)
    {
      current_condition = condition_comp( xxi, xet, xze );
      if ( current_condition > condition ) { condition = current_condition; }
    }

    if(metrics_request_flag & V_HEX_ODDY)
    {
      current_oddy = oddy_comp( xxi, xet, xze );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
    }
    

    // J(0,0,0)
    current_jacobian = edges[0] % (-edges[3] * edges[8]);
    if ( current_jacobian < jacobian )
      jacobian = current_jacobian;

    if(metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ))
    {
      det_sum += current_jacobian;
    }
    
    if(metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ))
    {
      current_scaled_jacobian = current_jacobian / 
        sqrt(length_squared[0] * length_squared[3] * length_squared[8]);
      if(current_scaled_jacobian < scaled_jacobian)
        shear = scaled_jacobian = current_scaled_jacobian;
    }

    if(metrics_request_flag & V_HEX_CONDITION)
    {
      current_condition = condition_comp( edges[0], -edges[3], edges[8] );
      if ( current_condition > condition ) { condition = current_condition; }
    }

    if(metrics_request_flag & V_HEX_ODDY)
    {
      current_oddy = oddy_comp( edges[0], -edges[3], edges[8] );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
    }

    if(metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ))
    {
      if(current_jacobian > VERDICT_DBL_MIN)
        current_shape = 3 * pow( current_jacobian, two_thirds) / 
          (length_squared[0] + length_squared[3] + length_squared[8]);
      else
        current_shape = 0;

      if(current_shape < shape) { shape = current_shape; }

    }

    // J(1,0,0)
    current_jacobian = edges[1] % (-edges[0] * edges[9]);
    if ( current_jacobian < jacobian )
      jacobian = current_jacobian;
    
    if(metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ))
    {
      det_sum += current_jacobian;
    }
    
    if(metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ))
    {
      current_scaled_jacobian = current_jacobian / 
        sqrt(length_squared[1] * length_squared[0] * length_squared[9]);
      if(current_scaled_jacobian < scaled_jacobian)
        shear = scaled_jacobian = current_scaled_jacobian;
    }

    if(metrics_request_flag & V_HEX_CONDITION)
    {
      current_condition = condition_comp( edges[1], -edges[0], edges[9] );
      if ( current_condition > condition ) { condition = current_condition; }
    }

    if(metrics_request_flag & V_HEX_ODDY)
    {
      current_oddy = oddy_comp( edges[1], -edges[0], edges[9] );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
    }
    
    if(metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ))
    {
      if(current_jacobian > VERDICT_DBL_MIN)
        current_shape = 3 * pow( current_jacobian, two_thirds) / 
          (length_squared[1] + length_squared[0] + length_squared[9]);
      else
        current_shape = 0;

      if(current_shape < shape) { shape = current_shape; }

    }

    // J(1,1,0)
    current_jacobian = edges[2] % (-edges[1] * edges[10]);
    if ( current_jacobian < jacobian )
      jacobian = current_jacobian;
    
    if(metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ))
    {
      det_sum += current_jacobian;
    }
    
    if(metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ))
    {
      current_scaled_jacobian = current_jacobian / 
        sqrt(length_squared[2] * length_squared[1] * length_squared[10]);
      if(current_scaled_jacobian < scaled_jacobian)
        shear = scaled_jacobian = current_scaled_jacobian;
    }
    
    if(metrics_request_flag & V_HEX_CONDITION)
    {
      current_condition = condition_comp( edges[2], -edges[1], edges[10] );
      if ( current_condition > condition ) { condition = current_condition; }
    }
    
    if(metrics_request_flag & V_HEX_ODDY)
    {
      current_oddy = oddy_comp( edges[2], -edges[1], edges[10] );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
    }

    if(metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ))
    {
      if(current_jacobian > VERDICT_DBL_MIN)
        current_shape = 3 * pow( current_jacobian, two_thirds) / 
          (length_squared[2] + length_squared[1] + length_squared[10]);
      else
        current_shape = 0;

      if(current_shape < shape) { shape = current_shape; }

    }


    // J(0,1,0)
    current_jacobian = edges[3] % (-edges[2] * edges[11]);
    if ( current_jacobian < jacobian )
      jacobian = current_jacobian;
    
    if(metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ))
    {
      det_sum += current_jacobian;
    }
    
    if(metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ))
    {
      current_scaled_jacobian = current_jacobian / 
        sqrt(length_squared[3] * length_squared[2] * length_squared[11]);
      if(current_scaled_jacobian < scaled_jacobian)
        shear = scaled_jacobian = current_scaled_jacobian;
    }

    if(metrics_request_flag & V_HEX_CONDITION)
    {
      current_condition = condition_comp( edges[3], -edges[2], edges[11] );
      if ( current_condition > condition ) { condition = current_condition; }
    }
    
    if(metrics_request_flag & V_HEX_ODDY)
    {
      current_oddy = oddy_comp( edges[3], -edges[2], edges[11] );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
    }

    if(metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ))
    {
      if(current_jacobian > VERDICT_DBL_MIN)
        current_shape = 3 * pow( current_jacobian, two_thirds) / 
          (length_squared[3] + length_squared[2] + length_squared[11]);
      else
        current_shape = 0;

      if(current_shape < shape) { shape = current_shape; }

    }

    // J(0,0,1)
    current_jacobian = edges[4] % (-edges[8] * -edges[7]);
    if ( current_jacobian < jacobian )
      jacobian = current_jacobian;
    
    if(metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ))
    {
      det_sum += current_jacobian;
    }
    
    if(metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ))
    {
      current_scaled_jacobian = current_jacobian / 
        sqrt(length_squared[4] * length_squared[8] * length_squared[7]);
      if(current_scaled_jacobian < scaled_jacobian)
        shear = scaled_jacobian = current_scaled_jacobian;
    }

    if(metrics_request_flag & V_HEX_CONDITION)
    {
      current_condition = condition_comp( edges[4], -edges[8], -edges[7] );
      if ( current_condition > condition ) { condition = current_condition; }
    }
    
    if(metrics_request_flag & V_HEX_ODDY)
    {
      current_oddy = oddy_comp( edges[4], -edges[8], -edges[7] );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
    }

    if(metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ))
    {
      if(current_jacobian > VERDICT_DBL_MIN)
        current_shape = 3 * pow( current_jacobian, two_thirds) / 
          (length_squared[4] + length_squared[8] + length_squared[7]);
      else
        current_shape = 0;

      if(current_shape < shape) { shape = current_shape; }

    }

    // J(1,0,1)
    current_jacobian = -edges[4] % (edges[5] * -edges[9]);
    if ( current_jacobian < jacobian )
      jacobian = current_jacobian;
    
    if(metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ))
    {
      det_sum += current_jacobian;
    }
    
    if(metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ))
    {
      current_scaled_jacobian = current_jacobian / 
        sqrt(length_squared[4] * length_squared[5] * length_squared[9]);
      if(current_scaled_jacobian < scaled_jacobian)
        shear = scaled_jacobian = current_scaled_jacobian;
    }
    
    if(metrics_request_flag & V_HEX_CONDITION)
    {
      current_condition = condition_comp( -edges[4], edges[5], -edges[9] );
      if ( current_condition > condition ) { condition = current_condition; }
    }
    
    if(metrics_request_flag & V_HEX_ODDY)
    {
      current_oddy = oddy_comp( -edges[4], edges[5], -edges[9] );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
    }

    if(metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ))
    {
      if(current_jacobian > VERDICT_DBL_MIN)
        current_shape = 3 * pow( current_jacobian, two_thirds) / 
          (length_squared[4] + length_squared[5] + length_squared[9]);
      else
        current_shape = 0;

      if(current_shape < shape) { shape = current_shape; }

    }

    // J(1,1,1)
    current_jacobian = -edges[5] % (edges[6] * -edges[10]);
    if ( current_jacobian < jacobian )
      jacobian = current_jacobian;
    
    if(metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ))
    {
      det_sum += current_jacobian;
    }
    
    if(metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ))
    {
      current_scaled_jacobian = current_jacobian / 
        sqrt(length_squared[5] * length_squared[6] * length_squared[10]);
      if(current_scaled_jacobian < scaled_jacobian)
        shear = scaled_jacobian = current_scaled_jacobian;
    }
    
    if(metrics_request_flag & V_HEX_CONDITION)
    {
      current_condition = condition_comp( -edges[5], edges[6], -edges[10] );
      if ( current_condition > condition ) { condition = current_condition; }
    }
    
    if(metrics_request_flag & V_HEX_ODDY)
    {
      current_oddy = oddy_comp( -edges[5], edges[6], -edges[10] );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
    }

    if(metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ))
    {
      if(current_jacobian > VERDICT_DBL_MIN)
        current_shape = 3 * pow( current_jacobian, two_thirds) / 
          (length_squared[5] + length_squared[6] + length_squared[10]);
      else
        current_shape = 0;

      if(current_shape < shape) { shape = current_shape; }

    }

    // J(0,1,1)
    current_jacobian = -edges[6] % (edges[7] * -edges[11]);
    if ( current_jacobian < jacobian )
      jacobian = current_jacobian;
    
    if(metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ))
    {
      det_sum += current_jacobian;
    }
    
    if(metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ))
    {
      current_scaled_jacobian = current_jacobian / 
        sqrt(length_squared[6] * length_squared[7] * length_squared[11]);
      if(current_scaled_jacobian < scaled_jacobian)
        shear = scaled_jacobian = current_scaled_jacobian;
    }
    
    if(metrics_request_flag & V_HEX_CONDITION)
    {
      current_condition = condition_comp( -edges[6], edges[7], -edges[11] );
      if ( current_condition > condition ) { condition = current_condition; }
    }
    
    if(metrics_request_flag & V_HEX_ODDY)
    {
      current_oddy = oddy_comp( -edges[6], edges[7], -edges[11] );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
    }

    if(metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ))
    {
      if(current_jacobian > VERDICT_DBL_MIN)
        current_shape = 3 * pow( current_jacobian, two_thirds) / 
          (length_squared[6] + length_squared[7] + length_squared[11]);
      else
        current_shape = 0;

      if(current_shape < shape) { shape = current_shape; }

    }

    if(metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ))
    {
      if(det_sum > VERDICT_DBL_MIN && rel_size_error != VERDICT_TRUE)
      {
        double tau = det_sum / ( 8 * detw );
        metric_vals->relative_size_squared = (VERDICT_REAL)std::min( tau*tau, 1.0/tau/tau);
      }
      else
        metric_vals->relative_size_squared = 0.0;
    }

    // set values from above calculations
    if(metrics_request_flag & V_HEX_JACOBIAN)
      metric_vals->jacobian = (VERDICT_REAL)jacobian;

    if(metrics_request_flag & V_HEX_SCALED_JACOBIAN)
      metric_vals->scaled_jacobian = (VERDICT_REAL)scaled_jacobian;

    if(metrics_request_flag & V_HEX_CONDITION)
      metric_vals->condition = (VERDICT_REAL)(condition/3.0);

    if(metrics_request_flag & V_HEX_SHEAR)
    {
      if( shear < VERDICT_DBL_MIN )  //shear has range 0 to +1
        shear = 0;
      metric_vals->shear = (VERDICT_REAL)shear;
    }

    if(metrics_request_flag & V_HEX_SHAPE)
      metric_vals->shape = (VERDICT_REAL)shape; 
  
    if(metrics_request_flag & V_HEX_SHAPE_AND_SIZE)
      metric_vals->shape_and_size = (VERDICT_REAL)(shape * metric_vals->relative_size_squared);
  
    if(metrics_request_flag & V_HEX_SHEAR_AND_SIZE)
      metric_vals->shear_and_size = (VERDICT_REAL)(shear * metric_vals->relative_size_squared);
  
    if(metrics_request_flag & V_HEX_ODDY)
      metric_vals->oddy = (VERDICT_REAL)oddy;

    if(metrics_request_flag & V_HEX_STRETCH)
    {
      static const double HEX_STRETCH_SCALE_FACTOR = sqrt(3.0);
      double min_edge=length_squared[0];
      for(int j=1; j<12; j++)
        min_edge = std::min(min_edge, length_squared[j]);

      double max_diag = diag_length(1, coordinates);
        
      metric_vals->stretch = (VERDICT_REAL)(HEX_STRETCH_SCALE_FACTOR * ( safe_ratio( sqrt(min_edge), max_diag) ));
    }
  }


  if(metrics_request_flag & V_HEX_DIAGONAL_RATIO)
    metric_vals->diagonal_ratio = v_hex_diagonal_ratio(num_nodes, coordinates);
  if(metrics_request_flag & V_HEX_MAX_DIAGONAL)
    metric_vals->max_diagonal = v_hex_max_diagonal(num_nodes, coordinates);
  if(metrics_request_flag & V_HEX_MIN_DIAGONAL)
    metric_vals->min_diagonal = v_hex_min_diagonal(num_nodes, coordinates);
  if(metrics_request_flag & V_HEX_DIMENSION)
    metric_vals->dimension = v_hex_dimension(num_nodes, coordinates);
  if(metrics_request_flag & V_HEX_DISTORTION)
                metric_vals->distortion = v_hex_distortion(num_nodes, coordinates);

  
  //take care of any overflow problems
  //aspect
  if( metric_vals->aspect > 0 )
    metric_vals->aspect = (VERDICT_REAL) std::min( metric_vals->aspect, VERDICT_DBL_MAX );
  else
    metric_vals->aspect = (VERDICT_REAL) std::max( metric_vals->aspect, -VERDICT_DBL_MAX );

  //skew
  if( metric_vals->skew > 0 )
    metric_vals->skew = (VERDICT_REAL) std::min( metric_vals->skew, VERDICT_DBL_MAX );
  else
    metric_vals->skew = (VERDICT_REAL) std::max( metric_vals->skew, -VERDICT_DBL_MAX );

  //taper
  if( metric_vals->taper > 0 )
    metric_vals->taper = (VERDICT_REAL) std::min( metric_vals->taper, VERDICT_DBL_MAX );
  else
    metric_vals->taper = (VERDICT_REAL) std::max( metric_vals->taper, -VERDICT_DBL_MAX );

  //volume
  if( metric_vals->volume > 0 )
    metric_vals->volume = (VERDICT_REAL) std::min( metric_vals->volume, VERDICT_DBL_MAX );
  else
    metric_vals->volume = (VERDICT_REAL) std::max( metric_vals->volume, -VERDICT_DBL_MAX );

  //stretch
  if( metric_vals->stretch > 0 )
    metric_vals->stretch = (VERDICT_REAL) std::min( metric_vals->stretch, VERDICT_DBL_MAX );
  else
    metric_vals->stretch = (VERDICT_REAL) std::max( metric_vals->stretch, -VERDICT_DBL_MAX );

  //diagonal ratio
  if( metric_vals->diagonal_ratio > 0 )
    metric_vals->diagonal_ratio = (VERDICT_REAL) std::min( metric_vals->diagonal_ratio, VERDICT_DBL_MAX );
  else
    metric_vals->diagonal_ratio = (VERDICT_REAL) std::max( metric_vals->diagonal_ratio, -VERDICT_DBL_MAX );

  //dimension
  if( metric_vals->dimension > 0 )
    metric_vals->dimension = (VERDICT_REAL) std::min( metric_vals->dimension, VERDICT_DBL_MAX );
  else
    metric_vals->dimension = (VERDICT_REAL) std::max( metric_vals->dimension, -VERDICT_DBL_MAX );

  //oddy
  if( metric_vals->oddy > 0 )
    metric_vals->oddy = (VERDICT_REAL) std::min( metric_vals->oddy, VERDICT_DBL_MAX );
  else
    metric_vals->oddy = (VERDICT_REAL) std::max( metric_vals->oddy, -VERDICT_DBL_MAX );

  //condition
  if( metric_vals->condition > 0 )
    metric_vals->condition = (VERDICT_REAL) std::min( metric_vals->condition, VERDICT_DBL_MAX );
  else
    metric_vals->condition = (VERDICT_REAL) std::max( metric_vals->condition, -VERDICT_DBL_MAX );

  //jacobian
  if( metric_vals->jacobian > 0 )
    metric_vals->jacobian = (VERDICT_REAL) std::min( metric_vals->jacobian, VERDICT_DBL_MAX );
  else
    metric_vals->jacobian = (VERDICT_REAL) std::max( metric_vals->jacobian, -VERDICT_DBL_MAX );

  //scaled_jacobian
  if( metric_vals->scaled_jacobian > 0 )
    metric_vals->scaled_jacobian = (VERDICT_REAL) std::min( metric_vals->scaled_jacobian, VERDICT_DBL_MAX );
  else
    metric_vals->scaled_jacobian = (VERDICT_REAL) std::max( metric_vals->scaled_jacobian, -VERDICT_DBL_MAX );

  //shear
  if( metric_vals->shear > 0 )
    metric_vals->shear = (VERDICT_REAL) std::min( metric_vals->shear, VERDICT_DBL_MAX );
  else
    metric_vals->shear = (VERDICT_REAL) std::max( metric_vals->shear, -VERDICT_DBL_MAX );

  //shape
  if( metric_vals->shape > 0 )
    metric_vals->shape = (VERDICT_REAL) std::min( metric_vals->shape, VERDICT_DBL_MAX );
  else
    metric_vals->shape = (VERDICT_REAL) std::max( metric_vals->shape, -VERDICT_DBL_MAX );

  //relative_size_squared
  if( metric_vals->relative_size_squared > 0 )
    metric_vals->relative_size_squared = (VERDICT_REAL) std::min( metric_vals->relative_size_squared, VERDICT_DBL_MAX );
  else
    metric_vals->relative_size_squared = (VERDICT_REAL) std::max( metric_vals->relative_size_squared, -VERDICT_DBL_MAX );

  //shape_and_size
  if( metric_vals->shape_and_size > 0 )
    metric_vals->shape_and_size = (VERDICT_REAL) std::min( metric_vals->shape_and_size, VERDICT_DBL_MAX );
  else
    metric_vals->shape_and_size = (VERDICT_REAL) std::max( metric_vals->shape_and_size, -VERDICT_DBL_MAX );

  //shear_and_size
  if( metric_vals->shear_and_size > 0 ) metric_vals->shear_and_size = (VERDICT_REAL) std::min( metric_vals->shear_and_size, VERDICT_DBL_MAX );
  else
    metric_vals->shear_and_size = (VERDICT_REAL) std::max( metric_vals->shear_and_size, -VERDICT_DBL_MAX );
  
  //distortion
  if( metric_vals->distortion > 0 )
    metric_vals->distortion = (VERDICT_REAL) std::min( metric_vals->distortion, VERDICT_DBL_MAX );
  else
    metric_vals->distortion = (VERDICT_REAL) std::max( metric_vals->distortion, -VERDICT_DBL_MAX );

}



